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Abstract 

The  manipulation  of  two  dimensional  elevation  maps  is  an  important  part  of  digital  car¬ 
tography.  In  many  situations,  these  maps  are  computed  by  interpolating  sparse  data  such  as 
isolated  elevation  points  obtained  from  stereo  matching.  In  this  paper,  we  present  a  surface 
interpolation  algorithm  based  on  variational  splines  which  is  well  suited  to  massively  parallel 
computers.  Using  multiresolution  parallel  relaxation,  we  can  efficiently  compute  the  interpo¬ 
lated  surface  and  also  have  local  control  over  its  continuity  and  smoothness.  We  apply  this 
technique  to  sparse  elevation  data  and  to  elevation  contours,  and  show  how  to  add  realistic 
fractal  detail  through  stochastic  relaxation.  We  also  present  a  multiresolution  decomposition 
algorithm  and  a  fast  parallel  3-D  rendering  algorithm. 

1  Introduction 

Digital  cartography  is  an  emerging  field  which  concerns  itself  with  the  application  of  com¬ 
puters  to  the  creation  and  manipulation  of  maps  [Hanson88].  The  data  being  manipulated 
can  be  represented  in  a  variety  of  formats;  of  these,  the  digital  terrain  model  (DTM)  is  one 
of  the  most  common  and  most  general  {USGS87].  A  DTM  represents  the  data  using  a  two- 
dimensional  array  which  encodes  some  local  feature  such  as  elevation,  terrain  classification, 
or  cultural  status. 
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The  creation  of  digital  terrain  models  often  involves  the  interpolation  of  sparse  elevation 
data  obtained  through  field  surveys  or  aerial  photogrammetry  [SlarnaSO].  Other  common 
operations  include  the  alignment  of  data  taken  from  different  samples  or  sources,  and  the 
automatic  correction  of  errors  or  biases  in  the  data.  A  variety  of  techniques  have  been 
developed  for  doing  surface  interpolation  [Fuchs77,  Bartels87,  Cendes87].  In  this  paper,  we 
will  introduce  a  new  class  of  interpolators  based  on  variational  splines  [Ahlberg67].  These 
interpolators  are  extremely  flexible  and  allow  local  control  over  quantities  such  as  continuity, 
orientation,  and  smoothness.  The  variational  splines  are  formulated  using  an  optimization 
framework  which  specifies  the  desirable  properties  for  the  interpolated  solution  (e.g.,  the 
closeness  of  fit.  to  the  data  and/or  the  smoothness  of  the  solution). 

The  variational  splines  that  result  from  this  framework  can  be  computed  in  two  differ¬ 
ent  ways.  The  first  involves  triangulating  the  two-dimensional  domain  of  the  digital  terrain 
model  using  the  sparse  data  points  as  corners  and  solving  for  the  coefiicient  of  the  polynomial 
(spline)  patches  on  each  triangle.  The  second  involves  using  a  fine  rectangular  discretization 
which  corresponds  to  the  final  resolution  of  the  map.  This  approach  allows  us  to  combine 
information  from  various  sources  (e.g.,  low-resolution  elevation  maps  obtained  through  pho¬ 
togrammetry  with  high-precision  survey  measurements)  and  to  have  local  control  over  the 
model  properties.  This  latter  approach  is  the  one  which  we  adopt  in  this  paper. 

To  convert  the  analytic  (continuous)  equations  describing  the  variational  spline  into  a 
discrete  set  of  equations,  we  use  either  finite  difference  or  finite  element  analysis.  The  result¬ 
ing  set  of  equations  can  then  be  solved  using  parallel  relaxation.  These  iterative  algorithms 
solve  the  system  of  equations  by  repeatedly  setting  each  grid  point  to  a  weighted  average  of 
its  neighbors.  They  are  thus  extremely  well  suited  to  massively  parallel  architectures  such 
as  the  Connection  Machine'^^  [Hillis85]. 

Unfortunately,  the  convergence  of  these  relaxation  algorithms  can  be  quite  slow.  In 
this  paper,  we  show  how  multiresolution  methods  such  as  multigrid  and  preconditioning 
with  hierarchical  basis  functions  can  be  used  to  accelerate  this  convergence.  The  resulting 
algorithms  have  a  performance  which  makes  them  suitable  for  real  applications. 

In  addition  to  developing  an  eflBcient  surface  interpolation  algorithm,  we  examine  other 
parallel  algorithms  for  map  manipulation.  We  show  how  a  stochastic  (random)  generalization 
of  surface  interpolation  can  be  used  to  add  fractal  detail  to  our  surface.  This  same  approach 
can  also  be  used  to  quantify  the  uncertainty  in  our  map  estimates.  We  examine  how  to 
provide  a  multiresolution  description  of  the  terrain  model,  which  may  be  useful  in  the  process 
of  extracting  a  reduced  description  based  on  prominent  features.  We  also  develop  a  fast 
parallel  rendering  algorithm  for  displaying  three-dimensional  views  of  the  surface. 

The  algorithms  described  in  this  paper  have  all  been  implemented  on  the  Connection 
Machine.  Their  suitability,  however,  is  not  limited  to  just  fine-grained  (SIMD)  architectures 
[PotterSS].  Because  of  their  regular  structure,  our  algorithms  are  well  suited  to  medium¬ 
grained  (MIMD)  architectures  and  even  vector  or  array  processors.  Some  of  these  algorithms 
have  also  been  implemented  on  standard  RISC  architectures  with  good  performance.  Devel- 
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oping  massively  parallel  algorithms  thus  allows  us  to  use  a  wide  range  of  current  architectures 
while  anticipating  the  increasing  parallelism  of  future  computing  hardware. 


1.1  Overview 

This  paper  starts  in  Section  2  with  a  review  of  surface  interpolation  using  variational  splines. 
We  present  a  mathematical  formulation  of  the  problem  and  discuss  surface  interpolation  vs. 
surface  approximation.  We  also  show  how  to  achieve  local  control  over  both  smoothness  and 
continuity.  Our  variational  formulation  is  compared  to  a  Bayesian  (probabilistic)  formula¬ 
tion.  In  Section  3,  we  derive  the  discrete  equations  for  the  surface  (elevation  map)  and  show 
how  they  can  be  solved  using  local  relaxation.  We  discuss  two  different  acceleration  methods: 
multigrid  relaxation,  and  conjugate  gradient  descent  preconditioned  with  hierarchical  basis 
functions.  In  Section  4,  we  discuss  how  interpolation  can  be  used  for  amplification,  i.e.,  to 
increase  the  resolution  of  the  data  and  to  fill  in  unknown  arecis.  We  also  show  how  to  add 
fractal  detail  to  the  surface  to  increase  the  realism  and  how  to  use  such  random  samples  to 
obtain  a  Monte  Carlo  esimate  of  the  local  uncertainty  (variance)  in  the  solution.  Section  5 
examines  the  opposite  problem,  that  of  reducing  a  terrain  map  to  a  more  abstract  or  com¬ 
pact  form.  We  discuss  how  a  relative  multiresolution  representation  can  be  used  to  aid  this 
processes.  Section  6  describes  an  efficient  parallel  algorithm  for  rendering  3-D  views  of  the 
terrain  model  with  optional  texture  or  intensity  mapping  from  aerial  images.  We  conclude 
the  paper  with  a  discussion  of  the  advantages  and  disadvantages  of  the  two-dimensional 
terrain  model  representation. 


2  Surface  interpolation  with  variational  splines 

The  problem  of  interpolating  a  two-dimensional  surface  such  as  an  elevation  map  through  a 
set  of  data  points  is  usually  miderconstrained,  i.e.,  there  are  many  possible  surfaces  which 
pass  through  the  given  points.  One  way  of  resolving  this  problem  is  to  heuristically  choose 
a  good  interpolating  algorithm,  e.g.,  to  triangulate  the  domain  and  use  piecewise  linear 
patches.  Another  possibility  is  to  cast  surface  interpolation  as  an  optimization  problem, 
e.g.,  to  maximize  the  smoothness  of  the  surface  while  minimizing  the  error  of  the  fit  to  the 
data  points.  This  latter  approach,  which  is  called  variational  spline  fitting  [Ahlberg67],  is 
the  one  we  adopt  in  this  paper. 

To  formally  characterize  the  optimization  problem,  we  use  theory  [Tikhonov??]. 

This  mathematical  technique  imposes  a  weak  smoothness  constraint  on  the  possible  solu¬ 
tions.  The  functional  to  be  minimized 

£(/)  =  £d(/;  (Pi))  +  A£.(/)  (1) 

is  thus  the  weighted  sum  of  two  terms:  the  data  compatibility  constraint  £d{f‘,  {Pi})  and  the 
smoothness  constraint  £s{f)-  The  regularization  parameter  A  is  used  to  adjust  the  closeness 
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of  the  fit  between  the  surface  and  the  data.  Regularization  has  recently  been  applied  to  a 
wide  range  of  problems  in  low-level  computer  vision  [PoggioSSa,  TerzopoulosSGa]. 

2.1  Interpolation  vs.  approximation 

The  data  compatibility  constraint  measures  the  distance  between  the  collection  of  depth 
values  {p,-}  =  {(u,-,  u,-,  d;)}  and  the  interpolated  surface  f{u,v).  In  its  simplest  form,  this 
constraint  is  a  weighted  squares  energy  measure 

^d(/;  {pi})  = 

where  the  weights  c;  are  inversely  related  to  the  variance  of  the  measurements  (c,-  = 

This  weighting  of  the  data  allows  us  to  combine  both  low-resolution  data,  such  as  that 
obtained  from  large-scale  photogrammetry,  with  high-resolution  data  obtained  from  field 
surveys  or  more  detailed  photosurveys.  It  also  allows  us  to  take  into  account  the  inherently 
uncertain  nature  of  elevation  measurements. 

If  we  wish  the  surface  to  pass  exactly  through  the  data  points  (i.e.,  to  have  an  inter¬ 
polating  instead  of  an  approximating  solution),  we  can  set  the  weights  Ci  to  a  very  large 
value,  or  equivalently,  we  can  set  the  regularization  parameter  A  — +  0.  This  does  not  pose 
any  significant  complications  in  our  framework  (it  is  easy  to  represent  sufficiently  small  or 
large  values  using  floating  point  numbers),  but  it  may  slow  down  the  convergence  of  certain 
algorithms  such  as  conjugate  gradient  if  we  do  not  use  the  correct  preconditioning  (Section 
3).  In  practice,  we  almost  never  desire  a  true  interpolated  solution  since  our  data  points  are 
never  without  error.* 

Choosing  an  optimal  value  of  A  is  an  active  area  of  research  in  regularization  theory.  A 
powerful  method  for  choosing  this  parameter  is  called  generalized  cross-validation  [Craven79]. 
It  involves  minimizing  the  distance  between  each  data  point  d,-  and  the  spline  approximation 
fit  to  the  remaining  data  points.  This  approach  is  quite  expensive,  since  for  each  sample 
value  of  A,  n  spline  fits  must  be  calculated,  where  n  is  the  number  of  data  points.  A 
simpler  approach  is  simply  to  adjust  A  until  the  variance  of  the  residuals  r,-  =  /(u, •,?;,•)  —  d,- 
matches  that  of  the  noise  model  (erf).  A  related  method  is  process  whitening  [Maybeck82, 
Marroquin85],  which  chooses  a  A  that  makes  the  residuals  as  uncorrelated  as  possible.  A 
number  of  other  approaches  to  this  problem,  including  a  new  method  based  on  maximum 
likelihood  estimation,  are  discussed  in  [Szeliski88]. 

*We  could  also  formulate  the  regularization  problem  as  min/  Saif)  that  |/(u;,t)))  —  d.j  <  e  for  all 
:  [PoggioSSa].  However,  problems  of  this  sort  are  more  difficult  to  solve  and  the  utility  of  this  formulation 
seems  questionable. 
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2.2  Smoothness  and  continuity 


The  smoothness  constraint  (or  stabilizer^  in  regularization  theory)  is  a  functional  or  norm 
of  f{u,v)  that  encodes  the  variation  in  the  surface.  Two  examples  of  possible  smoothness 
functionals  are  the  membrane  model 

W)  =  ^j  /  (/u  +  fv)  dv, 

which  is  a  small  deflection  approximation  of  the  surface  area,  and  the  thin  plate  model 

&(/)  =  5  /  /  (/;„ + VI. + fi)  du  dv, 

which  is  a  small  deflection  approximation  of  the  surface  curvature  (the  subscripts  indicate 
partial  derivatives)  [TerzopoulosSGb].  These  two  models  can  be  combined  into  a  single  func¬ 
tional 

^s(/)  =  ^// -  'r(u,v)][/„^  +  /v]  +  tCu,  v)[/i  +  +  fl]}  dudv  (3) 

where  p(u,v)  is  a  “rigidity”  function,  and  t(u,v)  is  a  “tension”  function  .  The  rigidity  and 
tension  functions  can  be  used  to  allow  depth  (p(u,v)  =  0)  and  orientation  (t(u,v)  =  0) 
discontinuities.  The  above  smoothness  constraint  corresponds  to  a  “generalized  piecewise 
continuous  spline  under  tension”  [TerzopoulosSGb]. 

To  see  the  effects  of  using  various  smoothness  constraints,  consider  the  nine  data  points 
shown  in  Figure  la.  The  interpolated  solution  using  a  membrane  model  is  shown  in  Figure 
lb.  Note  how  the  surface  has  visible  peaks  and  dips  corresponding  to  the  location  of  the  data 
points.  In  practice,  the  membrane  interpolator  is  not  sufficiently  smooth  for  cartographic 
applications.  The  interpolated  solution  using  a  thin  plate  model  is  shown  in  Figure  Ic.  This 
model  may  sometimes  prove  to  be  too  smooth,  but  is  often  a  good  choice.  The  controlled- 
continuity  thin-plate  is  shown  in  Figure  Id.  Note  that  a  depth  discontinuity  has  been 
introduced  along  the  left  edge  and  an  orientation  discontinuity  along  the  right. 


2.3  Local  control 


The  stabilizer  £’s(/)  described  by  (3)  is  an  example  of  the  more  general  controlled-continuity 
constraint 


W)  =  I 

^  m=0  i+fc=m 


m\ 


j\k\ 


d”^f(u,v) 

duWv^- 


2 

du  dv 


(4) 


This  formulation  allows  us  to  have  local  control  over  the  continuity  and  the  smoothness  of 
the  interpolated  solution.  The  values  for  the  continuity  and  smoothness  can  be  derived  from 
other  cartographic  features  such  as  rivers,  peaks  or  cultural  features.  We  will  present  some 
graphical  examples  of  this  local  control  in  Section  4. 
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Figure  1:  Sample  data  points  and  interpolated  solutions 
(a)  sample  data  points  (b)  membrane  interpolant  (c)  thin  plate  interpolant  (d)  thin  plate 
with  discontinuities 
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Discontinuities  in  the  elevation  map  (tears)  or  in  its  orientation  (creases)  can  be  achieved 
by  locally  setting  twi  =  rti2  =  0  or  =  0,  respectively.  These  discontinuities  usually  occur 
along  lines  or  curves  such  as  cliffs  (for  tears)  or  ravines  and  ridges  (for  creases).^  In  the 
discrete  formulation  of  our  problem  (Appendix  A),  these  features  are  encoded  using  binary 
line  and  crease  variables. 

Smooth  variations  in  the  correspond  to  local  variation  in  surface  roughness.^  These 
can  be  obtained  from  a  knowledge  of  the  local  topography  (mountainous,  rolling,  flat,  etc.). 
We  can  also  derive  these  parameters  through  a  statistical  analysis  of  our  terrain  model 
[Anderssen74,  Pentland84,  Hewett86].  In  the  discrete  formulation,  local  smoothness  control 
is  achieved  by  making  the  w^'s  discrete  fields  coincident  with  the  elevation  map  (Appendix 
A).  The  controlled-continuity  constraint  (4)  can  also  be  extended  to  model  anisotropic  sur¬ 
faces  [Boult86],  i.e.,  surfaces  such  as  parallel  mountain  ranges  whose  correlation  or  structure 
is  not  rotationally  symmetric.  In  our  framework,  this  could  be  achieved  by  using  a  separate 
coefficient  in  front  of  each  partial  derivative. 

Lastly,  we  can  incorporate  non-linear  constraints  into  our  regularization  framework.  For 
example,  we  could  constrain  the  elevation  map  to  have  a  local  maximum  at  a  known  peak 
location.  We  could  also  constrain  certain  areas  of  the  map  (such  as  man-made  areas  or 
lakes)  to  be  flat.  These  constraints  result  in  discrete  equations  that  are  no  longer  linear 
and  are  not  directly  amenable  to  the  same  optimization  techniques  which  we  study  in  this 
paper.  However,  since  the  algorithms  which  we  develop  are  iterative,  the  extension  to  non¬ 
linear  constraints  is  not  too  difficult  [Platt88].  This  should  be  a  fruitful  direction  for  future 
research. 


2.4  Relationship  to  Bayesian  estimation 


The  energy-based  framework  for  surface  interpolation  can  be  related  to  a  probabilistic 
(Bayesian)  framework  through  the  use  of  Gibbs  distributions  [Geman84].  In  the  Bayesian 
framework,  surface  interpolation  is  cast  as  an  optimal  estimation  problem.  A  prior  models 
provides  a  statistical  description  of  the  surface  which  we  are  trying  to  reconstruct. 
A  sensor  model,  P{{pi}\f),  describes  the  noisy  process  by  which  the  data  points  were  sam¬ 
pled  from  the  surface.  From  these  two  distribution,  we  can  compute  the  posterior  model, 
P(/|{p,-}),  using  Bayes’  Rule 


P{f\{pi}) 


p{{pmp{f) 


(5) 


where  P({p,})  is  used  to  normalize  the  posterior  distribution.  We  can  then  compute  the 
Maximum  a  Posteriori  (MAP)  estimate  by  finding  the  value  of  /  which  maximizes  (5).'' 

^The  problem  of  automatically  detecting  discontinuities  in  surfaces  has  been  extensively  studied  in  the 
field  of  computer  vision  [Marroquin84,  Blake87,  TerzopouIos86b,  Leclerc89]. 

^This  is  equivalent  to  having  a  spatially  varying  A. 

■^Other  estimates  can  also  be  obtained  from  the  posterior  distribution  [Marroquin85,  Szeliski88],  but  for 
our  application  they  are  equivalent  to  the  MAP  estimate. 
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We  can  relate  the  Bayesian  framework  to  the  energy-based  regularization  framework  by 


setting 

P{f)  ccexp(-54/)) 

(6) 

and 

^({P.-}I/)  «  exp  (-^d(/;  {p.})) . 

(7) 

For  this  choice  of  prior  and 

sensor  models,  we  find  that 

-P(/|{Pi})  «  exp  (-£(/)) 

(8) 

and  that  maximizing  the  posterior  probability  ■P(/|{pi})  is  equivalent  to  minimizing  the 
energy  S{f).  The  MAP  estimate  therefore  corresponds  to  the  usual  minimum  energy  solution 
obtained  from  regularization. 

This  observation  has  been  made  before  in  the  regularization  literature  [KimeldorfTO],  but 
it  has  not  been  exploited  to  a  great  degree.  One  way  we  can  make  use  of  this  relationship 
between  energies  and  log  probabilities  is  to  build  better  sensor  models.  For  example,  the 
probability  distribution  (7)  corresponding  to  our  weighted  squares  data  constraint  (2)  is 


i  v27r(j,' 


[di  - 

2af  )  • 


(9) 


From  this,  we  see  that  our  data  constraint  corresponds  to  an  assumption  that  each  elevation 
di  is  obtained  by  sampling  the  surface  /  at  location  (u,-,  u,)  and  adding  independent  Gaussian 
noise  with  variance  aj.  More  sophisticated  sensor  models  which  model  gross  errors  or  three- 
dimensional  uncertainty  in  position  can  be  developed  and  integrated  into  the  variational 
framework  by  setting  the  data  constraint  energy  to  the  negative  log  of  the  sensor  probability 
distribution  [SzeliskiSS] . 

The  prior  model  obtained  by  exponentiating  a  stabilizer  of  the  form  (4)  is  a  correlated 
Gaussian  field  whose  spectrum  is  determined  by  the  order  of  the  stabilizer  (i.e,,  the  choice  of 
Wm^s)  [Szeliski87].  For  the  membrane  and  thin  plate  models,  the  prior  model  is  a  stochastic 
fractal,  i.e.,  a  random  surface  which  is  self-affine  over  scale  [Mandelbrot82,  Voss85].  While 
this  observation  may  sound  counter-intuitive,  it  provides  a  statistical  model  for  the  surfaces 
which  we  are  interpolating.  Having  such  a  model  allows  us  to  generate  typical  samples  from 
this  distribution  which  we  can  use  to  check  the  validity  of  our  prior  model  [Geman84]. 

Under  the  Bayesian  interpretation,  the  posterior  distribution  P{f\{pi})  describes  a  family 
of  random  surfaces  drawn  from  the  prior  model  which  are  consistent  with  the  observed 
data.  The  most  likely  sample  from  this  distribution,  the  MAP  estimate,  corresponds  to 
the  minimum  energy  solution  of  (1),  as  shown  in  Figure  Id.  A  typical  sample  from  this 
distribution,  which  can  be  generated  using  stochastic  relaxation  (Section  3),  is  shown  in 
Figure  2.  We  can  use  such  random  samples  to  add  realism  to  our  interpolated  surfaces  and 
to  calculate  the  uncertainty  in  our  surface  estimates  (Section  4). 
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Figure  2;  Random  sample  from  posterior  distribution 

3  Parallel  relaxation  algorithms 

3.1  Finite  element  discretization 

To  minimize  (1),  we  apply  the  finite  element  method,  which  provides  a  systematic  approach 
to  the  discretization  and  solution  of  variational  spline  problems  [TerzopoulosS3].  We  dis¬ 
cretize  f{u,v)  on  a  regular  fine-grained  mesh  of  nodal  variables  which  coincides  with  the 
digital  terrain  model  we  wish  to  produce.  The  advantage  of  such  a  representation  is  that  it 
is  independent  of  the  location  of  the  data  points.  We  can  thus  add  more  points  or  combine 
data  from  different  resolutions  without  changing  the  algorithm.  It  is  also  easier  to  represent 
discontinuities  and  local  variations  in  smoothness  in  this  form.  The  regular  fine-grained  na¬ 
ture  of  the  mesh  leads  naturally  to  simple  local  parallel  algorithms  which  can  execute  on  a 
massively  parallel  array  of  processors.® 

Using  a  triangular  conforming  element  for  the  membrane  and  a  non- conforming  rectan¬ 
gular  element  for  the  thin  plate  [Terzopoulos83],  we  can  derive  the  energy  equations 

-E's(x)  =  —  -h  —  X{j)  ]  (10) 

(«M) 

^The  alternatives  to  our  fine-grained  discretization  include  finite  elements  defined  over  a  data-dependent 
triangulation  of  the  domain  [Fuchs77,  Cendes87]  and  kernel  splines  [BoultSGj.  Non-regular  (adaptive)  grids 
could  also  be  used  to  give  increased  resolution  where  needed. 
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(11) 


for  the  membrane  (the  subscripts  indicate  spatial  position)  and 
■£'s(x)  =  —h 

{'<:)  +  2{x{^ij+\  —  Xij+i  —  Xi+i,j  + 

+  (x.',;+i  -  2xij  + 

for  the  thin  plate,  where  h  =  |Au|  =  |Ayj  is  the  size  of  the  mesh  (isotropic  in  u  and  u). 
These  equations  hold  at  the  interior  of  the  surface.  Near  border  points  or  discontinuities 
some  of  the  energy  terms  are  dropped  or  replaced  by  lower  continuity  terms  (see  Appendix 
A).  The  equation  for  the  data  compatibility  energy  is  simply 

^d(x,  d)  =  i  Ci,j{xij  -  dij)^  (12) 

with  Cij  =  0  at  points  where  there  is  no  input  data. 

If  we  concatenate  all  the  nodal  variables  into  one  vector  x,  we  can  write  the  prior 

energy  model  as  one  quadratic  form 

=  ^x^AsX  (13) 

where  x  is  the  vector  of  nodal  variables,  i.e.,  x  =  {f{hi,  hj)}.  The  stiffness  matrix  As  is 
extremely  sparse,  having  at  most  13  entries  per  row.  Similarly,  the  data  constraint  in  (2) 
can  be  written  as 

£^d(x,  d)  =  i(x  -  d)^Ad(x  -  d),  (14) 

where  d  is  a  zero-padded  vector  of  elevation  values,  and  the  diagonal  matrix  Aj  has  entries 
C{  where  data  points  coincide  with  nodal  variables  and  zeros  elsewhere. 

Using  (13)  and  (14),  we  write  the  combined  energy  (1)  in  discrete  form, as 

E{x)  =  ix^Ax  —  x^b  +  k  (15) 

where  is  a  constant, 

A  =  A  As  -f  Ad,  and  b  =  Add. 

This  energy  function  has  a  minimum  at  x  =  x*,  the  solution  to  the  linear  system  of  algebraic 
equations 

Ax  =  b.  (16) 

3.2  Deterministic  and  stochastic  relaxation 

In  principle,  the  solution  to  (16)  can  be  found  using  either  direct  or  iterative  numerical 
methods.  Direct  methods  are  impractical  for  solving  large  systems  eissociated  with  fine 
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meshes  because  of  excessive  storage  requirements.®  Iterative  methods,  on  the  other  hand, 
do  not  require  any  additional  storage  and  are  amenable  to  parallel  implementations.  Here, 
we  will  study  four  iterative  methods:  Gauss-Seidel  relaxation,  weighted  Jacobi  relaxation, 
conjugate  gradient  descent,  and  stochastic  relaxation. 

Gauss-Seidel  relaxation  is  in  principle  a  sequential  algorithm  where  the  nodal  variables 
Xi  are  updated  one  at  a  time.  The  new  value  for  x,-  is  chosen  to  minimize  the  local  energy 

E{xi)  =  +  I  Y,  1  (1"^) 

2  J 

where  a,j  are  the  entries  of  A,  and  A;  expresses  the  fact  that  the  a,j  are  nonzero  only  for 
certain  neighbors  of  node  i  (see  Appendix  A  for  details).  When  using  Gauss-Seidel  relaxation, 
we  choose  the  new  node  value  which  minimizes  this  local  energy 


Xf  = 


bj  y  j  GijXj 

K  / 


(18) 


(for  pure  interpolation,  we  set  xf  =  di  at  points  coincident  with  data).  We  can  thus  rewrite 
the  local  energy  as 

£l(xO  =  ^a.,(xi-.Tt)2  +  ^;.  (19) 

Gauss-Seidel  relaxation  is  well  suited  to  sequential  architectures,  since  it  converges  faster 
than  the  parallel  Jabobi  algorithm  (see  below).  On  a  parallel  architecture,  we  can  update 
nodes  simultaneously  as  long  as  no  two  nodes  are  neighbors  of  each  other  (j  ^  iV,).  For  the 
membrane  interpolator,  this  can  be  achieved  by  using  a  red-black  coloring  (updating  the  red 
squares  on  a  checkerboard  in  parallel,  then  the  black  squares)  [BriggsS?].  For  the  thin  plate, 
we  have  to  use  at  least  5  colors,  and  it  is  often  easier  to  code  the  algorithm  using  9  colors 
arranged  in  a  3  x  3  array.  On  a  machine  like  the  Connection  Machine,  this  can  result  in  an 
algorithm  where  each  Gauss-Seidel  cycle  is  9  times  slower  than  a  Jacobi  cycle. 

Weighted  Jacobi  relaxation  is  the  fully  parallel  version  of  Gauss-Seidel.  To  ensure  that 
the  iterative  algorithm  converges,  we  take  a  smaller  step  than  we  would  in  Gauss-Seidel 


xf  -  Xi  +  (b{-  Y 

\  ) 


(20) 


where  the  neighborhood  N\  now  contains  the  point  i.  The  underrelaxation  parameter  w 
is  chosen  so  that  the  algorithm  convergences.  A  good  choice  for  u>  can  be  computed  by 
examining  the  eigenvalues  of  the  iteration  matrix  [YoungTl,  BriggsST].  The  parallel  version 
of  (20)  can  be  written  using  matrix  notation, 

xfc+i  ^  x'-'  +  u;D-^(b  -  Ax'^),  (21) 

®For  an  N  X  N  image,  we  require  O(N^)  storage  and  0(A'‘')  operations  for  a  direct  solution. 


11 


where  D  is  the  matrix  of  diagonal  elements  from  A  (the  superscripts  indicate  the  iteration 
number).  The  term  r*'  =  b  —  Ax*'  is  called  the  residual  and  measures  the  current  difference 
between  the  right  and  left  hand  sides  of  (16). 

While  both  of  these  algorithms  do  a  good  job  when  the  data  is  reasonably  dense,  their 
performance  degrades  severely  as  the  distance  between  data  points  increases  (see  [Szeliski89a] 
for  examples).  A  faster  algorithm  is  conjugate  gradient  descent  [PressSG].  In  this  algorithm, 
we  choose  a  descent  direction  p  and  take  an  optimally  sized  step  in  this  direction 

x*'+*  =x*'  +  Q*'p*'.  (22) 

The  step  size  a*  is  chosen  to  minimize  the  global  energy  (15),  and  can  be  shown  to  be 

a*  =  p*'^r*'/p*'^Ap*'. 

Initially,  the  direction  p  is  taken  to  be  the  residual  r,  which  is  also  the  negative  gradient  of 
the  energy  (15).  To  increase  the  convergence  of  the  algorithm,  we  require  that  successive 
directions  be  conjugate,  i.e.,  that  p^^Ap*”*  =  0.  We  can  ensure  this  by  setting 

p*'  =  r*'-/3V-'  (23) 

with 

/3^  =  r*'^Ap*'-Vp*'"*^Ap^-*. 

As  we  can  see  from  these  operations,  this  algorithm  is  well  suited  to  massively  parallel  archi¬ 
tectures  such  as  the  Connection  Machine  which  can  both  compute  weighted  local  averages 
(e.g.,  Ax)  and  global  summations  (e.g.,  p^r).  A  full  description  of  the  conjugate  gradient 
algorithm  is  given  Table  1  near  the  end  of  this  section. 

The  three  algorithms  which  we  have  just  described  can  be  used  to  find  the  surface  which 
minimizes  the  energy  (15).  In  the  Bayesian  interpretation,  this  corresponds  to  the  MAP 
estimate  for  the  surface.  If  we  wish  instead  to  generate  a  typical  sample  from  the  posterior 
distribution,  we  can  use  the  Gibbs  Sampler  [Geman84].  This  algorithm  is  a  stochastic 
generalization  of  Gauss-Seidel  relaxation.  Updating  the  points  one  at  a  time,  we  choose  the 
new  value  for  z,-  from  the  local  Boltzmann  distribution  corresponding  to  (19) 

P{xi)  (X  exp  {-E{xi)/T)  a  exp  .  (24) 

This  distribution  is  a  Gaussian  with  mean  equal  to  the  deterministic  update  value  xf  and 
a  variance  equal  to  Tfau.  Thus,  the  Gibbs  Sampler  is  equivalent  to  the  usual  Gauss-Seidel 
relaxation  algorithm  with  the  addition  of  some  locally  controlled  Gaussian  noise  at  each 
step.’^  The  amount  of  roughness  can  be  controlled  with  the  “temperature”  parameter  T. 
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3.3  Multigrid  relaxation 


Although  the  above  iterative  algorithms  will  eventually  compute  the  desired  solution, 
the  convergence  may  be  unacceptably  slow  in  practice  (see  [Szeliski89a]  for  examples).  To 
accelerate  convergence,  we  can  use  coarse- to-fine  relaxation  on  a  multiresolution  pyramid 
(Figure  3).  The  problem  is  first  solved  on  a  coarser  (smaller)  mesh,  and  the  solution  is  used 
as  an  initial  condition  for  the  next  finer  level. 

This  simple  idea  is  often  adequate  to  get  a  rough  solution.  However,  to  achieve  maximum 
accuracy,  a  more  sophisticated  approach  called  multigrid  is  usually  required.  Multigrid  re¬ 
laxation  algorithms  [Hackbusch85,  Briggs87]  are  based  on  the  observation  that  local  iterative 
methods  are  good  at  reducing  the  high  frequency  components  of  the  interpolation  error,  but 
are  poor  at  reducing  the  low  frequency  components.  By  solving  a  related  problem  on  the 
coarser  grid,  this  low  frequency  error  can  be  reduced  more  quickly. 

To  develop  a  multigrid  algorithm,  several  components  must  be  specified:  the  method  used 
to  derive  the  energy  equations  at  the  coarser  levels  from  the  fine  level  equations;  a  restriction 
operation  that  maps  a  solution  at  a  fine  level  to  a  coarser  grid;  a  prolongation  operation  that 
maps  from  the  coarse  to  the  fine  level;  and  a  coordination  scheme  that  specifies  the  number 
of  iterations  at  each  level  and  the  sequence  of  prolongations  and  restrictions. 

A  sophisticated  version  of  multigrid  relaxation  was  developed  by  [Terzopoulos83].  In  his 
implementation,  the  coarser  meshes  are  coincident  with  the  original  mesh,  i.e.,  the  coarser 
level  nodes  are  a  subset  of  the  finer  level  nodes.  Simple  injection  (subsampling)  is  used 
for  restriction,  and  third-degree  Lagrange  interpolation  is  used  for  prolongation.  A  Full 
Multigrid  (FM)  method  is  used  to  coordinate  the  inter-level  interactions.  This  ensures  that 
the  coarse  levels  have  the  same  accuracy  as  the  fine  level  solutions. 

Our  own  experiments  with  multigrid  on  the  Connection  Machine  were  somewhat  dis¬ 
appointing.  The  convergence  rates  were  not  significantly  better  than  those  obtained  with 
simple  coarse- to-fine  relaxation.  This  problem  may  be  due  to  the  irregular  nature  of  the 
data  constraints  (multigrid  techniques  were  originally  developed  for  the  solution  of  uniform 
smooth  partial  differential  equations).  It  is  quite  possible  that  the  use  of  different  (locally 
adjusted)  interpolators  could  improve  the  convergence,  but  we  have  not  investigated  this 
possibility.  Instead,  we  have  been  developing  a  multiresolution  version  of  conjugate  gradient 
descent,  which  we  describe  next. 

3.4  Preconditioned  conjugate  gradient  descent 

The  conjugate  gradient  descent  algorithm  that  we  described  earlier  can  be  accelerated  dra¬ 
matically  using  the  hierarchical  bcisis  functions  developed  recently  by  [Yserentant86].  In 
this  approach,  the  usual  nodal  basis  set  x  is  replaced  by  a  hierarchical  basis  set  y.  Certain 

parallel  version  of  stochastic  relaxation  based  on  global  diffusion  is  also  possible  but  does  not  converge 
as  quickly  to  the  equilibrium  distribution  [Szeliski89b]. 
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elements  of  the  hierarchical  baisis  set  have  larger  support  than  the  nodal  basis  elements,  and 
this  allows  the  relaxation  algorithm  to  converge  more  quickly  when  using  the  hierarchical 
set. 

To  convert  from  the  hierarchical  to  the  nodal  basis  set,  we  use  a  simple  linear  (matrix) 
transform 

x  =  Sy,  (25) 

where  S  =  S1S2  . . .  Sl-i,  and  L  is  the  number  of  levels  in  the  hierarchical  basis  set.  Each  of 
the  sparse  matrices  S/  interpolates  the  nodes  at  level  /  +  1  to  level  I  and  adds  in  the  nodes 
corresponding  to  the  new  level.  The  columns  of  S  give  the  values  of  the  hierarchical  basis 
functions  at  the  nodal  variable  locations. 

In  his  paper,  Yserentant  uses  recursive  subdivision  of  triangles  to  obtain  the  nodal  basis 
set.  The  corresponding  hierarchical  basis  then  consists  of  the  top-level  (coarse)  triangulariza- 
tion,  along  with  the  subtriangles  that  are  generated  each  time  a  larger  triangle  is  subdivided. 
Linear  interpolation  is  used  on  a  triangle  each  time  it  is  subdivided.  We  can  generalize  this 
notion  to  arbitrary  interpolants  defined  over  a  rectangular  grid  [Szeliski89a] .  Each  node  in 
the  hierarchical  basis  is  assigned  to  the  level  in  the  multiresolution  pyramid  where  it  first 
appears  (the  solid  dots  in  Figure  3).  This  is  similar  to  the  multiresolution  pyramid,  except 
that  each  level  is  only  partially  populated,  and  the  total  number  of  nodes  is  the  same  in  both 
the  nodal  and  hierarchical  basis  sets.  To  fully  define  the  hierarchical  basis  set,  we  select  an 
interpolation  function  that  defines  how  each  level  is  interpolated  to  the  next  finer  level  before 
the  new  node  values  are  added  in. 

The  resulting  algorithms  for  mapping  between  the  hierarchical  and  nodal  basis  sets  are 
simple  and  efficient.  We  use 

procedure  S 

for  /  =  L  —  1  down  to  1 
for  i  6  Ml 
for  j  e  M 

u(z)  =  u(i)  -I-  iy(f;  j')u(j) 

end  S 

to  convert  from  the  hierarchical  to  the  nodal  basis  set.  In  this  procedure,  which  goes  from 
the  coarsest  level  (/  =  L)  to  the  finest  (/  =  1),  each  node  is  assigned  to  one  of  the  level 
collections  Mi.  Each  node  also  has  a  number  of  neighboring  nodes  W)  on  the  next  coarser 
level  that  contribute  to  its  value  during  the  interpolation  process.  The  w{i‘  j)  are  the  weight¬ 
ing  functions  that  depend  on  the  particular  choice  of  interpolation  function  (these  are  the 
off-diagonal  terms  in  the  S;  matrices). 

We  also  use  the  adjoint  of  this  operation 
procedure 
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for  /  =  1  to  L  —  1 
for  i  G  Mi 
for  j  G  Mi 

u{j)  =  u{j)  +  w{i-,j)u{i) 

end 

in  our  modified  conjugate  gradient  descent  algorithm. 

These  mapping  algorithms  are  easy  to  code  (once  the  interpolation  functions  have  been 
precomputed)  and  require  very  few  computations  to  perform.  On  a  serial  machine,  these 
procedures  use  0(n)  operations  (multiplications  and  additions),  where  n  is  the  number  of 
nodes.  On  a  parallel  machine,  0{L)  parallel  steps  are  required,  where  L  «  log  n  is  the 
number  of  levels  in  the  pyramid.  This  is  the  same  number  of  steps  as  is  needed  to  perform 
the  global  summations  (inner  products)  used  in  the  conjugate  gradient  algorithm.  Note  that 
although  we  have  defined  the  hierarchical  basis  over  a  pyramid,  it  can  actually  be  represented 
in  the  same  space  as  the  usual  nodal  basis,  and  the  transformations  between  bases  can  be 
accomplished  in  place. 

The  hierarchical  basis  set  allows  us  to  minimize  exactly  the  same  energy  as  the  one  we 
derived  from  the  discretization  on  the  finest  grid.  Substituting  x  =  Sy  into  (15),  we  obtain 
the  new  energy  equation 

E{y)  =  ■  iy^(S^AS)y  -  y^(S^b)  +  c 

=  -^y^Ay  -  y^b  +  c  (26) 

where  the  identifies  the  hierarchical  basis  vectors  and  matrices.  The  advantage  of  minimiz¬ 
ing  this  new  equation  is  that  the  condition  number  of  the  matrix  A  is  much  smaller  than 
that  of  the  original  matrix  A  [Yserentant86],  implying  much  faster  convergence  for  iterative 
algorithms  such  as  conjugate  gradient. 

Unfortunately,  the  A  matrix  is  not  as  sparse  as  the  original  matrix  A,  so  that  a  direct 
minimization  of  (26)  is  impractical.  Instead,  we  use  the  recursive  mappings  S  and  in 
conjunction  with  the  original  matrix  A  and  vector  b  to  compute  the  required  residuals  and 
inner  products.  The  resulting  conjugate  gradient  descent  algorithm  is  identical  to  the  usual 
single-level  algorithm  except  that  we  use  a  smoothed  version  of  the  residual  f  =  SS^r  to 
choose  the  new  direction  [Szeliski89a].  The  full  algorithm  is  shown  in  Table  1. 

The  above  algorithm  is  an  example  of  the  more  general  idea  of  preconditioning.  One 
common  application  of  this  idea  is  to  divide  the  residual  r  by  the  diagonal  of  the  A  matrix, 
i.e.,  to  use  r  =  D“^r.  This  is  equivalent  to  dividing  each  row  and  column  of  the  A  matrix  by 
the  square  root  of  the  diagonal  element  (A  =  D~*/^AD~^/^).  The  new  stiffness  matrix  A 
has  a  lower  condition  number,  especially  when  the  data  coupling  terms  are  very  stiff  (large 
c,-  or  small  A).  The  resulting  conjugate  gradient  algorithm  looks  a  lot  like  weighted  Jacobi 
with  the  added  conjugation  in  the  descent  direction.  We  can  combine  the  hierarchical  basis 
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0.  initialize  Tq  and  Fq  using  6.  —  8.  and  set  Pq  =  fq 

1.  wjt  =  Apjt 

2.  af  =  pk  ■  w/,  =  pjApjt 

3.  of  =  Pk-rk  =  plTk 

4. T  ak  =  a^lcc^ 

5.  xjt+i  =  xjt  +  ajtpjt 

6.  rjt+i  =  b-Axjt+i 

7.  Tjt+i  =  (preconditioning  step) 

8.  rjt+i  =  Srjt+i  (preconditioning  step) 

.  9.  =  Ffc+i-WA- 

104 

lb  pjt+i  =  — /Jjt+i  Pa- 

12.  loop  to  1. 

t  A£(x  +  op)  =  ^p^Ap  -  ap^r 
J  (f  -  /?P)^Ap  =  0 

Table  1:  Conjugate  gradient  descent  algorithm 

preconditioner  with  the  inverse  diagonal  preconditioner  (using  either  or  to 

obtain  an  algorithm  which  performs  well  even  when  the  data  coupling  is  stiff. 

3.5  Comparison  of  algorithms 

To  compare  the  performance  of  the  various  deterministic  relaxation  algorithm,  we  ran  these 
algorithms  on  the  data  set  shown  in  Figure  5.  First,  we  computed  the  optimal  solution  x" 
using  500  iterations  of  the  hierarchical  basis  conjugate  gradient  algorithm  {L  =  4).  Then, 
for  each  of  the  algorithms  tested,  we  computed  the  root  mean  squared  (RMS)  error 

tk  =  |xa  -x*l/Vn 

as  a  function  of  the  number  of  iterations  k. 

Figure  4  shows  the  relative  performance  of  the  relaxation  algorithms,  as  well  as  the  effects 
of  varying  the  number  of  smoothing  levels  in  hierarchical  basis  conjugate  gradient.  The 
topmost  curve  corresponds  to  weighted  Jacobi  relaxation  (cu  =  0.25),  which  is  the  slowest  of 
the  algorithms  tested.  The  second  curve  is  for  Gauss-Seidel,  which  initially  converges  faster 
but  then  also  tails  off.  Conjugate  gradient  does  better  and  continues  to  converge  towards 
the  solution  at  a  reasonable  rate.  Using  preconditioning  with  hierarchical  basis  functions 
dramatically  speeds  up  the  convergence  of  conjugate  gradient.  Initially,  the  speedup  is 
proportional  to  the  number  of  smoothing  levels.  As  the  iterations  progress,  however,  the 
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iterations 


Figure  4:  Convergence  plot  for  surface  interpolation 
The  top  three  curves  are  for  Jacobi,  Gauss-Seidel,  and  conjugate  gradient.  The  lower  five 
curves  are  for  hierarchical  basis  conjugate  gradient  with  different  numbers  of  smoothing  levels 
L  (note  that  HBCG  with  L  =  I  corresponds  to  conjugate  gradient).  The  fastest  convergence 
is  for  HBCG  with  Z  =  3  or  i  =  4. 
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versions  using  an  intermediate  number  of  smoothing  levels  (Zr  =  3  and  Z,  =  4)  have  the 
best  performance.  We  believe  that  the  optimum  number  of  levels  to  be  used  is  related  to 
the  density  of  the  underlying  data  points,  and  we  are  currently  investigating  methods  for 
automatically  adjusting  the  smoothing  based  on  the  local  structure  of  the  data. 


4  Data  amplification 

The  deterministic  relaxation  algorithms  which  we  have  developed  can  be  used  to  interpolate 
smooth  terrain  models  from  sparse  irregular  data.  Our  stochastic  relaxation  algorithm  can 
be  used  to  generate  random  terrain  models  which  fit  the  data.  We  can  use  both  of  these 
techniques  for  data  amplification,  i.e.,  to  produce  a  denser  (and  more  voluminous)  description 
of  our  terrain  from  a  sparse  set  of  control  points. 

4.1  Interpolation  from  sparse  data 

The  variational  splines  which  we  use  for  surface  interpolation  are  well  suited  to  interpolating 
sparse  irregular  elevation  data.  An  example  of  such  data  is  that  obtained  from  matching 
isolated  features  in  a  stereo  pair  of  aerial  images.  Figure  5a  shows  one  of  the  two  intensity 
images  which  was  input  to  the  STEREOSYS  stereo  matching  system  developed  at  SRI 
[HannahSS].  Figure  5b  shows  the  set  of  7696  matches  which  were  produced  by  this  algorithm 
in  a  sub-region  of  this  image.  These  points  were  then  given  as  sparse  elevation  constraints 
to  our  thin  plate  interpolating  program.  The  resulting  dense  elevation  map  is  shown  as  a 
contour  map  in  Figure  5c. 

Variational  splines  can  also  be  used  to  generate  a  full  digital  elevation  model  from  a  set  of 
contour  lines.  Figure  6a  shows  the  original  elevation  map  (from  a  1:250,000  scale  map  of  the 
Rockies  just  west  of  Denver)  digitized  on  a  256  x  256  grid.  In  this  figure,  brighter  intensities 
represent  higher  elevations,  and  the  contour  lines  are  drawn  at  100m  intervals.  This  map  Wcis 
subsampled  using  200m  contour  lines,  and  the  resulting  sparse  data  was  input  to  the  surface 
fitting  algorithm.  The  resulting  elevation  maps  are  shown  in  Figures  6b  (membrane),  6c  (thin 
plate),  and  6d  (mixed  membrane  /  thin  plate  model).  Although  it  is  difficult  to  see  from  the 
contour  maps,  the  membrane  interpolator  =  {OjljO})  is  not  smooth  enough 

(locally,  it  acts  like  a  “soap  film”),  while  the  thin  plate  interpolator  ({ti^o,  ,  ^^2}  =  {0, 0, 1}) 

is  too  smooth  (it  creates  a  large  “dip”  in  the  middle  of  the  map).  Using  a  blend  of  these  two 
models  =  {0,0.0063,1})  results  in  a  better  reconstruction. 

We  can  also  use  additional  (non-topographic)  information  to  control  and  enhance  the 
surface  interpolation  process.  An  example  of  such  information  is  the  hydrography  (streams, 
rivers,  and  other  drainage)  that  is  often  available  in  conjunction  with  digital  elevation  models 
[USGS86].  Figure  7a  shows  a  256  x  256  section  of  a  digital  elevation  model  of  the  Big  Bcisin 
area;  Figure  7b  shows  the  hydrography  for  the  same  region.  If  we  interpolate  a  surface 
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Figure  5:  Surface  interpolated  from  sparse  elevation  data 
(a)  sample  aerial  photograph  from  stereo  pair  (b)  matched  points  from  stereo  pair  (c)  inter¬ 
polated  thin  plate  solution  on  106  x  99  grid 
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directly  from  the  contours  shown  in  Figure  7c,  we  obtain  the  surface  shown  in  Figure  8a.  If, 
on  the  other  hand,  we  use  the  hydrography  to  introduce  creases  (discontinuities  in  the  first 
derivative),  we  obtain  the  more  realistic  (and  accurate)  map  shown  in  Figure  Sb.  A  similar 
effect  is  shown  in  Figures  8c  and  8d,  which  were  interpolated  from  the  2%  sampling  shown 
in  Figure  7d. 

For  interpolating  regularly  spaced  data  or  re-sampling  grids  for  rectification,  other  tech¬ 
niques  may  be  preferable  to  variational  splines.  Simple  bilinear  interpolation  or  cubic  in¬ 
terpolation  (Cendes87,  FoleyS7]  may  suffice.  In  any  event,  the  interpolation  of  regularly 
spaced  data  with  no  discontinuities  or  local  variation  can  always  be  expressed  as  a  linear 
shift-invariant  filtering  operation  [Oppe'nheim75,  Poggio85b].  This  results  in  a  simpler  and 
quicker  algorithm  than  the  one  we  use  in  this  paper.  Variational  splines  are  preferred, 
however,  when  there  are  significant  gaps  in  the  data  and  we  desire  local  control  over  the 
continuity. 

4.2  Adding  fractal  detail 

Stochastic  fractals  have  become  a  popular  technique  in  computer  graphics  for  adding  realistic 
detail  to  rendered  images  of  natural  terrain  [Fournier82,  Mandelbrot82].  Such  realism  is  use¬ 
ful  in  flight  simulators,  synthetic  scenes  used  in  advertising  and  movies,  and  as  backgrounds 
for  architectural  studies.  As  we  saw  in  Section  2.4,  a  random  sample  from  the  family  of 
surfaces  defined  by  the  variational  spline  energy  has  fractal  statistics  [Szeliski87].  In  Section 
3,  we  saw  how  a  simple  stochastic  extension  of  relaxation  (i.e.,  the  addition  of  controlled 
Gaussian  noise)  can  be  used  to  generate  fractal  surfaces  which  also  match  our  desired  data 
constraints. 

If  we  apply  our  stochastic  relaxation  algorithm  to  the  data  shown  in  Figure  6a,  we  obtain 
the  contour  map  shown  in  Figure  9.  Note  how  the  new  contour  lines  look  more  realistic  than 
the  smooth  ones  shown  in  Figure  6d.  A  shaded  rendering  of  this  same  surface  in  shown  in 
Figure  10.  We  can  also  use  our  stochastic  relaxation  with  additional  inter-level  consistency 
constraints  to  generate  a  “zoom  sequence”  which  reveals  more  detail  as  the  viewer  approaches 
the  surface  [Szeliski89b]. 


4.3  Illusion  vs.  reality 

The  surface  shown  in  Figures  9  and  10  is,  of  course,  just  one  possible  random  sample  from  a 
whole  family  of  surfaces  which  fit  the  data.  The  detail  in  the  surface  is  actually  “hallucinated” 
and  does  not  necessarily  correspond  to  the  true  detailed  shape  of  the  terrain.  For  certain 
applications  (fly-throughs,  mission  planning  and  simulation),  having  this  fictitious  detail 
may  be  advantageous.  For  others,  however,  this  detail  gives  a  false  sense  of  confidence  in 
the  terrain  shape  which  is  not  justified  by  the  available  data. 
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Figure  10:  Fractal  surface  (shaded  rendering) 
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From  a  Bayesian  perspective,  the  fractal  surface  of  Figure  9  is  a  typical  sample  from  the 
posterior  distribution,  whereas  Figure  6d  is  the  most  likely  sample  (and  also  the  mean  value). 
By  generating  a  number  of  random  samples,  we  can  actually  compute  the  local  variance  in 
the  posterior  estimate  at  each  point.®  These  confidence  measures  can  help  decide  where 
more  measurements  are  needed  or  how  much  to  trust  the  terrain  model  at  a  given  point. 

Displaying  the  variance  as  a  confidence  interval  around  each  surface  point  is  difficult.  In 
one  dimension,  we  can  draw  an  envelope  around  the  curve  showing  the  confidence  interval. 
In  two  dimensions,  we  could  use  color  coding  to  indicate  the  confidence.  If  our  computation 
and  rendering  hardware  were  fast  enough,  we  could  actually  display  many  random  samples 
in  succession.  The  resulting  surface  would  appear  to  “wiggle”  or  “shimmer”  in  areas  where 
it  is  most  uncertain.  Whether  this  additional  confidence  information  will  be  useful  for 
cartographic  applications  remains  to  be  seen. 


5  Data  reduction 

The  converse  of  the  data  amplification  process  discussed  in  the  previous  section  is  data 
reduction.  Here,  the  aim  is  to  go  from  a  complete  description  such  as  an  elevation  map  or 
terrain  model  to  a  more  parsimonious  description.  For  example,  we  may  wish  to  describe 
the  topography  by  the  location  of  simple  features  such  as  peaks,  ridges,  and  valleys.  Such 
sketch  maps  can  be  used  in  the  recognition  of  features  and  observer  location,  the  semantic 
interpretation  of  the  terrain  model,  and  the  quick  generation  of  rendered  views. 

The  general  problem  of  deriving  a  reduced  description  is  a  difficult  one  which  involves 
elements  of  both  signal  processing  and  knowledge-based  artificial  intelligence  [FischlerSQ].  In 
this  section,  we  will  sketch  out  a  more  modest  proposal  which  aims  to  decompose  the  terrain 
model  into  a  multiresolution  description.  This  work,  while  incomplete,  may  lead  to  a  novel 
technique  for  decomposing  a  surface  into  a  hierarchical  description. 


5.1  Relative  depth  representation 

The  multigrid  and  coarse-to-fine  algorithms  described  in  Section  3.3  keep  a  full  copy  of 
the  current  depth  estimate  at  each  level.  An  alternative  to  this  absolute  multiresolution 
representation  is  a  relative  representation  [SzeliskiSS].  In  this  representation,  each  level  yi 
encodes  details  pertinent  to  its  own  scale,  and  the  sum  of  all  the  interpolated  levels  provides 
the  overall  depth  map  estimate 

x  =  I]l/y/  (27) 

where  I;  is  the  interpolation  function  associated  with  each  level.  The  relative  representation 
is  thus  analogous  to  a  band-pass  image  pyramid  [Burt83,  Crowley82],  while  the  absolute 

®This  is  a  Monte  Carlo  approach  to  computing  a  distribution. 
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representation  is  similar  to  a  low-pass  pyramid.  In  the  relative  representation,  each  level  in 
the  pyramid  hats  its  own  eissociated  smoothness  energy 

E'Av')  =  (28) 

The  data  compatibility  is  meeisured  using  the  summed  depth  map  x  and  the  data  points  d 

£d(x,d)  =  i(x~dfAd(x-d).  (29) 

Using  y  =  . . .  y^^]^  to  denote  the  concatenation  of  the  individual  state  vectors,  I  = 

[I^ . .  .1'^]  for  the  concatenated  interpolation  matrices,  and 

■A'  0  •••  0 

0  A^  0 

As  = 

.0  0  •••  Af 

to  denote  the  composite  spline  model  matrix,  we  can  write  the  overall  energy  function  as 

E(y)  =  l;£;(y')  +  £d(^.d) 

1=1 

=  jy^'A.y  +  5(iy  -  d)"'Ai(iy  -  d) 

=  jy’'Ay-y’'b  +  c  (31) 

where 

A  =  As+FAdi  and  b  =  i^Aad.  (32) 

The  minimum  energy  solution  can  be  calculated  from  this  quadratic  form  by  solving 

Ay  =  b.  (33) 

A  number  of  iterative  relaxation  techniques  could  be  used  to  implement  the  minimiza¬ 
tion  of  the  composite  energy  (31).  For  example,  we  could  set  the  energy  gradient  of  (31) 
to  0  at  each  level  independently  using  the  current  value  of  x  to  compute  V£d(x,d),  and 
then  recompute  the  summed  representation  x.  Because  of  the  tight  coupling  between  the 
levels,  however,  this  approach  (Gauss-Seidel  or  Jacobi)  does  not  work  very  well.  Instead,  we 
have  found  that  using  conjugate  gradient  descent  performs  satisfactorily.  To  implement  this 
algorithm,  we  must  compute  the  current  residual 

r  =  Asy  +  FAd(iy  -  d). 

This  is  achieved  by  first  computing  x  =  ly,  then  subtracting  the  data  points  d  and  weighting 
the  error  by  Ad.  This  data-fit  residual  is  then  converted  into  the  relative  representation  by 


(34) 
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multiplying  by  (the  adjoint  of  the  interpolation  function),  and  the  smoothness  component 
of  the  residual  (A^y^)  is  then  added.  A  similar  approach  is  used  to  compute  the  sparse  matrix 
product  Ap,  where  p  is  the  current  direction  vector.  While  this  algorithm  does  not  converge 
as  fast  as  hierarchical  basis  conjugate  gradient,  it  does  produce  a  multiresolution  description 
while  solving  the  surface  interpolation  problem. 

To  reduce  the  amount  of  storage  required  for  the  interpolation  matrices  and  the  amount 
of  computation  performed  when  interpolating,  we  can  use.  a  recursive  structure  for  the  in- 
terpolants.  We  define 

z'  =  y'  +  lLi2'-'  (35) 

as  the  absolute  representation  at  level  /  (with  x  =  Zi),  where  l{_]  is  the  interpolation 
function  between  levels  /  —  1  and  /.  The  structure  of  our  relative  representation  algorithm 
now  looks  very  similar  to  that  of  the  hierarchical  basis  conjugate  gradient,  except  that  the 
smoothness  constraint  is  computed  at  each  level  separately,  and  each  level  in  the  pyramid  is 
fully  populated.  The  relative  representation  thus  has  more  state  than  the  original  fine-level 
system  of  equations. 

Designing  a  set  of  spline  energies  for  each  level  that  decomposes  x  into  a  reasonable 
multiresolution  description  while  maintaining  the  faithfulness  of  the  energy  function  (31)  to 
the  original  energy  (15)  is  a  nontrivial  task.  The  discussion  in  [SzeliskiSS]  gives  conditions 
for  the  approximation  to  hold.  It  is  often  easier  to  take  an  indirect  approach  to  designing 
these  energies  using  Fourier  analysis,  where  instead  of  defining  the  interpolation  problem  as 
a  functional  minimization  problem,  we  view  it  as  a  Bayesian  estimation  problem  and  try  to 
match  the  power  spectra  of  the  prior  models  [SzeliskiSS]. 

The  relative  representation  as  we  have  described  it  so  far  could  also  be  obtained  by 
filtering  the  fine-level  (composite)  representation.®  What  may  be  more  interesting  is  to 
equip  each  level  with  a  non-linear  smoothness  constraint,  e.g.,  one  that  favors  flat  or  zero 
surfaces.  Instead  of  making  the  cost  proportional  to  the  square  of  the  surface  magnitude, 
slope  or  curvature,  we  could  incur  a  fixed  cost  whenever  one  of  these  quantities  is  non-zero. 
For  example,  we  could  add  terms  of  the  form 

£.(x)  =  5;(1-«(j:,j))  (36) 

(*.i) 

This  would  favor  surfaces  which  are  zero,  and  higher  order  terms  could  be  used  to  favor 
piecewise  constant  or  piecewise  linear  solutions.^®  Combining  this  idea  with  a  multiresolution 
relative  representation,  we  would  obtain  an  interpolation  algorithm  which  distributes  the 
terrain  features  among  the  various  levels  in  a  non-linear  fashion.^*  While  this  is  still  far 
removed  from  a  true  semantic  interpretation  of  our  terrain,  it  may  provide  the  kind  of 
perceptual  grouping  which  is  a  necessary  precursor  of  object  identification  [Lowe85]. 

®This  is  because  the  levels  are  a  linear  function  of  the  input  data. 

^®This  idea  is  inspired  by  a  similar  constraint  in  Leclerc’s  work  on  minimal-length  encoding  [Leclerc89], 
where  the  optimal  polynomial  order  for  a  surface  patch  is  determined  by  penalizing  non-zero  coefficients. 

’■^By  contrast,  the  current  relative  representation  distributes  the  data  among  the  levels  in  a  “weak”  fashion, 
as  if  springs  of  varying  stiffness  were  attached  to  the  surface. 
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Figure  11:  Texture  mapped  rendering  of  Yosemite  valley 

6  Rendering 

Once  we  have  interpolated  our  cartographic  data  to  its  final  resolution,  we  may  wish  to 
display  this  information.  Cartographic  information  is  often  displayed  in  two-dimensional 
formats  such  as  contour  maps  (e.g.,  Figure  6a).  In  many  cases,  however,  a  three-dimensional 
perspective  view  may  be  more  useful  for  interpreting  or  analyzing  the  terrain. 

A  number  of  formats  can  be  used  for  three-dimensional  displays.  We  can  use  a  wireframe 
with  hidden  surface  removal  (Figure  1),  but  this  often  limits  the  available  resolution  of  the 
display.  We  can  also  use  a  simple  shading  model  to  generate  a  three-dimensional  view  (Figure 

10) .  Finally,  if  an  intensity  image  (aerial  photograph)  is  available  which  corresponds  to  the 
terrain  model,  we  can  use  texture  mapping  to  produce  a  realistic  view  of  the  scene  (Figure 

11) . 

The  computation  of  a  realistic,  properly  anti-aliased  rendering  can  be  quite  complicated 
and  time  consuming  [HansonSS].  For  many  applications,  however,  having  a  real-time  ren¬ 
dering  capability  greatly  enhances  the  interpretation  process.  This  is  especially  true  of 
the  real-time  rotation  of  three-dimensional  surfaces,  since  a  very  realistic  percept  of  shape 
can  be  obtained.^^  Traditionally,  such  animations  have  been  available  only  from  graphics 
workstations^^  or  flight  simulators.  In  our  own  research,  we  have  developed  a  rendering 

the  vision  literature,  this  is  called  the  kinetic  depth  effect. 

^^For  example,  the  top  of  the  line  IRIS  workstation  from  Silicon  Graphics  will  render  100,000  polygons 
per  second. 
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algorithm  for  the  Connection  Machine  which  gives  near  real-time  performance. 

6.1  A  parallel  rendering  algorithm 

The  key  to  developing  a  parallel  rendering  algorithm  for  the  Connection  Machine  is  to 
assign  one  processor  to  each  point  in  the  terrain  map  being  rendered.  Our  algorithm  also 
takes  advajitage  of  the  combining  routing  network  of  the  Connection  Machine  to  perform 
z-buffering.  A  detailed  description  of  the  algorithm  follows. 

The  input  to  our  rendering  algorithm  consists  of  two  images  (2-D  arrays):  a  floating  point 
image  of  elevation  values,  and  an  integer  image  of  intensities  to  be  texture  mapped.  These 
latter  intensities  can  be  true  intensities  obtained  from  aerial  photographs,  synthetic  intensi¬ 
ties  computed  using  a  shading  model,  or  a  simple  texture  such  as  a  grid.  The  elevation  map 
is  converted  into  screen  coordinates  (x,  y,  and  screen-based  depth  z)  using  a  user-supplied 
homogeneous  transformation  matrix  [Newman79]  (we  keep  an  extra  4  bits  of  precision  on 
the  screen  coordinates). 

Next,  we  optionally  blur  the  texture  (intensity)  image  to  reduce  the  aliasing  which  occurs 
when  multiple  points  map  to  the  same  screen  location.  To  do  this,  we  average  each  point 
with  its  nearest  neighbors  if  these  neighbors  map  to  the  same  or  nearby  points  on  the  screen 
(a  global  parameter  called  the  blur  radius  controls  the  nearness  criterion).  We  then  average 
the  points  from  the  previous  step  with  their  neighbors  two  pixels  away  if  they  are  still  within 
the  selected  screen  distance.  This  iteration  continues  until  the  alieising  has  been  eliminated. 

At  this  stage,  we  can  independently  draw  each  pixel  into  a  z-buffer  which  corresponds 
to  the  final  image  size.  To  do  this,  we  append  the  intensity  value  for  each  map  point  to 
its  depth  (z)  value,  i.e.,  we  form  a  long  (32-bit)  integer  which  contains  the  z  value  in  its 
high-order  bits  and  the  intensity  in  its  low-order  bits.  Using  the  global  routing  network  of 
the  Connection  Machine  with  max  combining^'*  we  can  directly  draw  each  pixel  into  the 
z-buffer  (pixels  falling  outside  the  z-buffer  are  automatically  clipped).  The  low-order  bits  of 
the  z-buffer  then  correspond  to  the  rendered  image. 

This  “simple”  rendering  is  fast,  but  it  often  leaves  gaps  in  the  image  where  a  map 
grid  square  covers  more  than  one  pixel.  While  this  fast  rendering  may  be  suitable  for 
quick  interactive  positioning  and  for  obtaining  a  depth  percept,  a  better  algorithm  must  be 
devised  for  final  renderings.  We  designed  our  “complex”  rendering  algorithm  by  assigning  a 
processor  to  each  square  interior  to  the  digital  map  (this  is  basically  equivalent  to  the  node 
per  processor  cissignment  with  some  replication  at  the  borders).  Each  processor  then  checks 
if  all  four  corners  of  the  square  map  to  the  same  screen  location.  If  so,  the  pixel  is  drawn;  if 
not,  the  square  is  subdivided  (refined)  before  drawing. 

The  subdivision  occurs  by  first  determining  the  maximum  number  of  screen  pixels  sepa¬ 
rating  the  square  corners  along  each  dimension 

”  m3,x(jx,-^.i —  x,-j]  -h  |yt+i,j  —  -}-  |y,q-ij+i  —  y,-j-|-i|,l) 

^^The  *Lisp  instruction  is  (+pset  :max  .  . .). 
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—  rnax(|a;,’j+i  "I"  |yt,j+i  yi'j'li  l^»+i,i+i  “t"  |yi+ij+i  yt+i,i|il)' 

We  then  subdivide  each  square  into  m,j  x  riij  pieces,  thus  guaranteeing  that  each  piece 
covers  at  most  one  pixel.  To  implement  the  subdivision,  each  square  processor  “grabs” 
,s,-j  =  mijTiij  processors  using  a  (scan!  [  Sij  '  +  !  !)  global  scan  to  determine  the  starting 
address  of  the  group  of  processors  it  will  use.  Often,  there  are  not  enough  processors  available 
in  the  rendering  processor  set,  in  which  case  the  complete  subdivision/rendering  stage  may 
take  several  iterations.  Each  group  of  processors  renders  its  assigned  square  in  parallel 
by  first  computing  the  screen  addresses  and  intensity  using  bilinear  interpolation  and  then 
calling  the  “simple”  rendering  procedure.’® 

The  performance  which  we  have  achieved  with  this  algorithm  is  quite  impressive.  On  our 
8K  Connection  Mzw;hine  running  unoptimized  *Lisp  code,  we  can  render  a  256  x  256  texture- 
mapped  image  in  0.45  seconds  for  the  simple  version  and  2.61  seconds  for  the  complex  fully 
interpolated  version  (these  figures  vary,  depending  on  the  particular  viewing  position). 

There  are  a  number  of  extensions  possible  to  our  rendering  algorithm  which  we  have 
not  implemented.  Instead  of  mapping  an  intensity  image  onto  the  terrain  map,  we  could 
map  the  surface  normals.  This  would  allow  us  to  Itiove  the  light  source  position  and  recom¬ 
pute  the  synthetic  shaded  image  in  real  time,  which  would  aid  in  the  perception  of  surface 
shape.’®  We  could  also  generalize  the  algorithm  to  render  arbitrary  3-D  surfaces  defined 
using  polygonal  (quadrilateral)  surface  patches.  This  would  only  require  a  different  front 
end  to  the  algorithm,  since  the  same  routines  could  be  used  for  coordinate  transformation, 
quadrilateral  subdivision,  and  z-bufFer  rendering.  Finally,  if  we  are  to  achieve  the  same  qual¬ 
ity  of  renderings  as  is  currently  available  with  the  SRI  Cartographic  Modeling  Environment 
[Hanson88],  we  would  have  to  add  a  multiresolution  pyramid  representation  of  the  terrain 
and  intensity  to  our  algorithm.  Integrating  this  representation  into  our  current  algorithm  is 
difficult,  since  it  is  based  on  projecting  each  map  point  into  the  rendered  image  rather  than 
ray  tracing  screen  points  back  to  the  map.  However,  a  recursive  rendering  algorithm  based 
on  the  idea  of  quadtrees  [Rosenfeld80]  may  be  feasible. 


7  Discussion 

The  approach  to  surface  interpolation  developed  in  this  paper  is  based  on  using  a  regular 
grid  representation  for  the  terrain  map.  This  representation  fits  in  well  with  the  usual 
desired  output  of  the  interpolation  stage.  It  is  limited,  however,  to  representing  single¬ 
valued  functions  such  as  terrain  (for  a  description  of  how  to  integrate  this  representation 
with  three-dimensional  models  of  buildings  and  objects,  see  [Hanson88,  BobickS9].). 

^^Anothet  alternative  which  we  considered  was  iteratively  splitting  each  square  into  two  until  it  is  small 
enough.  This  approach,  however,  requires  keeping  a  list  of  “active”  squares  which  have  not  yet  been  painted, 
and  ensuring  that  this  list  does  not  overflow  existing  memory  is  difficult. 

’®This  suggestion  was  made  by  Yvan  Leclerc. 
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The  spline  interpolator  defined  over  this  regular  grid  results  in  a  very  large  sparse  set  of 
equations  which  must  be  solved  iteratively.  Fortunately,  multiresolution  acceleration  tech¬ 
niques  (in  particular,  our  hierarchical  basis  conjugate  gradient  algorithm)  make  the  solution 
of  these  equations  feasible  in  a  reasonable  amount  of  time. 

These  iterative  relaxation  algorithms  are  inherently  parallel,  and  are  therefore  well  matched 
to  massively  parallel  architectures  such  as  the  Connection  Machine.  Multiresolution  algo¬ 
rithms,  however,  can  sometimes  underutilize  the  parallelism  in  such  an  architecture,  espe¬ 
cially  if  the  coarser  levels  of  the  pyramid  are  smaller  than  the  processor  array  size.^^  Because 
of  their  regularity,  our  parallel  algorithms  should  also  execute  efficiently  on  sequential,  vec¬ 
tor,  or  MIMD  (coarse-grained)  parallel  architectures.  The  key  to  the  widespread  availability 
and  acceptance  of  our  parallel  relaxation  algorithms  will  probably  lie  in  our  ability  to  auto¬ 
matically  covert  code  from  a  high-level  parallel  language  such  as  *Lisp  into  languages  such  as 
Common  Lisp,  C,  or  FORTRAN  than  can  be  compiled  on  more  conventional  architectures. 
This  promises  to  be  an  interesting  area  for  future  research. 


8  Conclusions 

In  this  paper,  we  have  demonstrated  how  variational  splines  can  be  used  for  performing 
surface  interpolation,  which  is  an  important  problem  in  digital  cartography.  Our  variational 
splines  are  defined  using  an  optimization  framework,  where  the  characteristics  of  the  desired 
solution  are  specified  rather  than  the  algorithm  used  to  compute  the  solution.  This  opti¬ 
mization  framework  is  useful  for  integrating  data  from  multiple  sources  and  resolutions,  and 
also  allows  us  to  locally  control  surface  characteristics  such  as  continuity  and  smoothness. 

To  compute  the  variational  spline,  we  use  a  fine,  regular  discretization  which  corresponds 
in  resolution  to  the  final  digital  map  we  wish  to  produce.  This  discretization  results  in  a 
simple  local  set  of  equations  and  is  well  suited  to  parallel  representations  and  algorithms.  To 
efficiently  solve  these  equations,  we  must  use  some  multiresolution  acceleration  technique.  In 
our  experience,  conjugate  gradient  descent  preconditioned  with  hierarchical  basis  functions 
has  proven  to  be  the  most  efficient  method. 

We  have  applied  our  surface  interpolation  algorithm  to  the  problem  of  data  amplification, 
i.e.,  obtaining  a  dense  cartographic  representation  from  a  sparse  set  of  samples.  Examples 
with  both  sparse  (punctate)  data  and  contour  lines  were  presented.  We  also  showed  how  to 
locally  control  the  surface,  for  example  by  introducing  creases  along  hydrographic  features. 
Natural  looking  fractal  detail  can  be  added  to  our  interpolated  solution  through  stochastic 
relaxation.  The  random  samples  thus  generated  can  be  used  to  quantify  the  local  uncertainty 
(variance)  in  our  terrain  model  estimates. 

The  problem  of  data  reduction  was  also  studied  briefly.  We  described  a  relative  repre¬ 
sentation  surface  interpolation  algorithm  which  decomposes  the  surface  into  a  hierarchy  of 

this  case,  the  combination  of  a  smaller  parallel  machine  and  virtual  processors  [TMC88]  may  actually 
be  a  boon. 
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scales,  and  suggested  how  this  could  be  extended  to  perform  a  more  non-linear  feature  ex¬ 
traction.  In  the  final  part  of  the  paper,  we  developed  a  parallel  three-dimensional  rendering 
algorithm  which  can  be  used  to  visualize  our  interpolated  surfaces. 

As  we  have  demonstrated  in  this  paper,  manipulating  terrain  maps  using  parallel  relax¬ 
ation  algorithms  gives  the  user  more  flexibility  in  integrating  data  and  adding  local  con¬ 
straints.  The  advent  of  parallel  architectures  and  the  development  of  fast  parallel  algorithms 
is  quickly  making  this  approach  practical  for  real  cartographies  applications. 
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A  Finite  element  equations 

In  this  Appendix,  we  present  the  discrete  implementation  of  the  energy  equations  used  in 
Section  3.1.  The  energy  equations  which  we  develop  are  formulated  in  terms  of  a  number 
of  discrete  fields  which  form  the  basic  data  structures  used  in  the  surface  interpolation 
algorithm: 


Xi,j 

surface  depth  (free  variables) 

dij 

depth  constraints 

depth  constraint  weights  (ct“^) 

U,j  > 

depth  discontinuities  [0,1] 

"ui 

orientation  discontinuities  [0,1] 

The  line  variables  Uj  and  m{j  are  located  on  a  dual  grid,  and  the  crease  variables  n,j  are 
coincident  with  the  depth  value  nodes. 

To  describe  the  discrete  version  of  the  smoothness  constraint,  we  define  a  number  of 
implicit  fields  which  are  not  actually  computed  by  the  algorithm  but  serve  only  as  a  notational 
shorthand.  First,  we  define  the  finite  differences 


Xt,j+2  j+l  XiJ. 

Next,  we  define  the  continuity  strengths 

K,i  =  (i-ki) 

=  (1 

^i,j  —  (1  “  ~  ~  j+i)(l  ~  j+i) 

These  continuity  strengths  determine  which  of  the  molecules  defined  by  the  . . . ,  x,""- 
finite  differences  are  active  (we  use  phantom  line  variables  around  the  edges  of  the  grid  to 
take  care  of  boundary  conditions). 
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Figure  12;  Continuity  strengths  and  computational  molecules 


The  interactions  between  the  depth  field,  the  line  and  crease  variables,  and  the  continuity 
strengths  are  shown  in  Figure  12.  Intuitively,  the  line  variables  disable  any  molecules  that 
are  cut  when  the  line  variable  is  turned  on.  The  crease  variables  disable  the  or 
molecules  whose  centers  are  coincident  with  the  crease,  and  disable  the  x^''-  molecules  if  two 
opposite  corners  are  creased. 

The  prior  energy  (weak  smoothness  constraint)  is  constructed  from  these  elements 


+  w,h-^  i)' 


(•-i) 


(37) 


where  h  =  Au  =  Av  is  the  grid  spacing.  The  weighting  functions  Wm,  which  come  from  the 
general  controlled-continuity  constraint  (4),  define  the  order  of  the  interpolator.  In  terms  of 
the  notation  used  by  Terzopoulos  [Terzopoulos86b]  which  is  shown  in  (3),  we  have 


{u;o,u;i,u;2}  =  {0,/)(u,u)[l  -  T(u,u)],/)(u,i;)r(n,i;)} 

(in  our  implementation,  we  use  the  line  and  create  variables  to  represent  the  discontinuities 
rather  than  spatially  varying  w-m)- 

The  equation  for  the  data  compatibility  constraint  is  the  same  as  we  presented  in  Section 
3.1,  namely 

^d(x,  d)  =  ^  (38) 

(«.i) 

with  Cij  =  being  the  inverse  variance  of  the  measurement  dij. 

The  stiffness  matrix  A  and  the  weighted  data  vector  b  can  be  obtained  from  the  discrete 
energy  equation 

F;(x)  =  £s(x)  +  £d(x,  d) 
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by  setting 


and. 


d'^Ejx) 

dx{jdxk,i  jj_Q 


b 


i,3 


dEjx) 

x=0 


(39) 

(40) 


The  A  matrix  and  b  vector  are  precomputed  before  the  relaxation  begins.  Since  A  has  at 
most  13  non-zero  entries  per  row,  we  can  store  these  coefficients  in  an  M  x  A  x  13  array, 
where  M  and  N  are  the  terrain  model  dimensions. 


Acknowledgements 

This  research  was  supported  by  a  Visiting  Scientist  position  at  the  Artificial  Intelligence 
Center  of  SRI  International.  I  would  like  to  thank  Martin  Fischler,  Yvan  Leclerc,  Lynn 
Quam  and  the  other  membeEof  the  Perception  Group  for  their  ideas  and  encouragement. 


39 


